Dynamics of F = 1 87 Rb condensates at finite temperatures 
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We investigate the dynamics of a F = 1 spinor Bose-Einstein condensate of 87 Rb atoms confined 
in a quasi-one-dimensional trap both at zero and at finite temperature. At zero temperature, we 
observe coherent oscillations between populations of the various spin components and formation of 
multiple domains in the condensate. We study also finite temperature effects in the spin dynamics 
taking into account the phase fluctuations in the Bogoliubov-de Gennes framework. At finite T, 
despite complex multi-domain formation in the condensate, population equipartition occurs. The 
length scale of these spin domains seems to be determined intrinsically by non-linear interactions. 

PACS numbers: 03.75.Mn, 03.75.Kk 



INTRODUCTION 

The seminal theory papers of T.-L. Ho [l| and T. Ohmi 
and K. Machida Q, as well as the experiments performed 
by the MIT group on optically trapped sodium Bosc- 
Einstein condensates (BEC) 3] have stimulated the de- 
velopment of a new interesting area of research in the 
field of multi-component quantum gases: BEC with spin 
degrees of freedom, i.e. the so called spinor condensates. 
The bosonic quantum field operator in such systems is no 
longer a scalar, but a vector, and depending on various 
parameters the system can be found in a coherent super- 
position or in an incoherent mixture of condensates with 
different spin components, that can exchange population 
depending on the intrinsic nonlinear coupling. The cou- 
plings between the different spin components, similarly 
as in the standard scalar BEC can be described by zero 
range potentials that however, acquire a matrix form and 
describe couplings that assure total spin conservation in 
elastic binary collisions |^y. 

It is worth mentioning that spinor condensates are 
closely related to the effective spin- 1/2 systems realized 
by radio-frequency coupling of the two hyperfine states 
m in such systems spin-waves above the critical 

temperature [|| and decoherence effects were observed 
0. From the theoretical point of view, collective modes 
in two-component BEC's have been studied in while 
an analysis on the formation of spin domains can be 
found e.g. in 

Spinor condensates are suitable systems to study var- 
ious quantum phenomena, that do not occur in single 
component BEC's. Equilibrium states of spinor conden- 
sates in an optical trap can exhibit magnetic ordering of 
various kinds. For instance, for sodium [F = 1), the free- 



dom of spin orientation leads to the formation of spin do- 
mains in an external magnetic field, which can be either 
miscible or immiscible with one another 0- Excitations 
in such systems may manifest the different spin char- 
acter; stability of topological excitations and textures, 
such as ordinary vortices, coreless vortices (Skyrmions) 
0,0], and t'Hooft-Polyakov monopoles 0], depends in 
a complex manner on the parameters of the system. The 
complexity of these systems becomes even more trans- 
parent in the limit of strongly correlated systems: ultra- 
cold Bose spin gases in optical lattices exhibit fascinating 
properties, including new types of quantum phases, such 
as a polar condensate phase, a condensate of singlet pairs, 
a crystal spin nematic phase, and a spin singlet crystal 
phase [r2l ]. 

The ground state, its magnetic properties and the low 
temperature thermodynamics of spinor condensates have 
been studied in several experiments: in 23 Na (F = 1) 
which has an anti-ferromagne tic grou nd state, in 87 Rb in 
the F = 1 spin state 15. 16. 17] which is ferromag- 
netic, as well as in the F = 2 spin state, which presents a 
rich ground state behavior [1J, [l7| • Recent experiments 
involving non-destructive imaging of magnetization of a 
spin F = 1 Bose gas of 87 Rb atoms [l8| open a new route 
to study magnetism of spinor condensates. Also, the re- 
cent experiment with chromium fl9j condensates opens 
the possibility to study even higher spin states (F = 3) 
with long range (dipolar) interactions and, presumably a 
much more complex phase diagram. 

In the last two years the focus of investigation has 
shifted towards two other very important aspects of the 
spinor BEC physics: the equilibrium properties at finite 
T, and dynamics of spinor BEC's. The thermodynamic 
properties of an ultra-cold spin Bose gas have been exper- 
imentally investigated in Ref. (2lJ. Recently, finite tern- 
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perature effects to describe the properties of the equilib- 
rium density distribution have been considered within the 
Hartree-Fock- Popov theory [2l| . Spin dynamics has been 
experimentally studied in 87 Rb condensates for F = 1 in 
[l5|. and more exhaustively for F = 2 in [3 EE E3- 
Coherent collisional spin dynamics in 87 Rb for F = 2 has 
been recently observed also in optical lattices - these 
experiments pave way towards the efficient creation of 
entangled atom pairs in optical lattices. 

Here we focus on the dynamics of F = 1 elongated 
spinor condensates at zero and finite temperatures. One 
of our major motivation to study the thermal effects on 
spin dynamics is to learn something about decoherence 
in multi-component systems. In the case F = 1 the in- 
ternal coupling of the spin components depends only on 
2 coefficients, while for F = 2 it depends on 3. Due to 
this simplicity, the case F = 1 allows for a better under- 
standing of the interplay between nonlinear interactions, 
spin couplings and temperature effects. We consider here 
the full coupled dynamical equations of the spin compo- 
nents, obtained within a mean- field framework, without 
any further approximation such as the Single Mode Ap- 
proximation (SMA) uty . or variational ansatz [23, 
that would mask some aspects of the complex dynami- 
cal evolution. We restrict our analysis to the quasi-lD 
case, when the condensates are kinematically frozen in 
the transverse directions. We investigate then the influ- 
ence of thermal effects on the spin dynamics. To this aim 
we use the Bogoliubov-de Gennes description of phonon 
modes and treat them as classical random fields, simi- 
larly as in Ref. j^. For quasi-lD condensates at finite 
T the main contribution to phonon fluctuations comes 
from the fluctuations of the condensate phase [27| , and 
this contribution is fully accounted in our simulations. 

The paper is organized as follows. First, in section II 
we introduce the model for T = 0. We describe some 
of the details of the numerical method in section III. In 
section IV we describe our results for T = 0. Section V 
we present a discussion of finite temperature effects. We 
conclude in section VI. 



DESCRIPTION OF THE SYSTEM 

The many-body Hamiltonian describing a F=l spinor 
condensate in absence of an external magnetic field is 
given by 

H = jd 3 r Ul (-^V 2 + vj\ * m (1) 

where v I / m (r) (3^) is the field operator that annihilates 
(creates) an atom in the 771-th hyperfine state \F = 1, m) 
at point r (m = 1,0,-1). The trapping potential 
Vext(r) is assumed harmonic and spin independent. The 



terms with coefficients Co and c 2 describe binary elas- 
tic collisions of spin-1 atoms in the combined symmet- 
ric channel of total spin and 2, and are expressed 
in terms of the s-wave scattering lengths ag and a 2 : 
c = 4nh 2 (a + 2a 2 )/3M and c 2 = 47r/7 2 (a 2 - a )/3M, 
where M is the atomic mass. F are the spin-1 matrices 
P, l2l| . The system is ferromagnetic if c 2 < ( 87 Rb), 
and anti- ferromagnetic if c 2 > ( 23 Na). 



l*o| 2 



The total number of atoms N = f dr(\*f?i\ 
|*_i| 2 ) and the total magnetization M = J dr(\^i\ 2 — 
l^-il 2 ) commute with the Hamiltonian, and are thus 
constants of motion. 

In the mean field approach, a condensate order param- 
eter is introduced for each magnetic sublevel ip m (r) = 
(^ m (r)), and by neglecting quantum fluctuations it 
yields the following energy functional 



d 3 r 



2M 



+ V'jV'm + yC^jFmfc • F j7 ^<0fc . (2) 

According to ihdip m /dt — 58/5^^, the coupled dy- 
namical equations for the spin components are obtained 

MM 



ih 



dt 



2M 



V 



eff 



ip m + c 2 T* 



(3) 



n m (r) = \ip m \ 2 is the density of the m-th component and 
77 = |?Ai| 2 + | V'o 1 2 + l^-i 1 2 is t ne total density normal- 
ized to the total number of atoms N. The population of 
the hyperfine state \l,m) is N m = J dr\ip m \ 2 such that 
N = Nt+No+N-t. We have defined T± ± = t/jfy* 



T * = 2tpiij}Qtf)-i, and the effective potentials that will 
determine the spatial dynamics of each component 



V±{' = V ex t + c ri + c 2 (±77i + 77 T n-x) 



eff _ 



V ext + C 71 + C 2 (l7l + ?7_l) . 



(4) 



Analogously to a spin-polarized condensate, the multi- 
component Gross-Pitaevskii equations J31 can be rewrit- 
ten in the form of continuity equations 



d 

dl 1 



Vj, 



5n m (r,t) . 



(5) 



where j m = h(ip^ n 'Vi/j m — ^ m 'Vip^ n )/2iM is the current, 
and Sn m (r, t) = —(2c%/K)Im(T m 'iJ) m ) is the rate of trans- 
fer of populations between spin components. Since the 
total number of atoms and magnetization are conserved, 
and the dynamical equations for m — ±1 are symmetric, 
it is verified that Sho = — 2Sh±i = (2c 2 /ih) (^pi^pQ 2 ip-i — 

In our calculation, we assume an axially-symmetric 
harmonic confinement V ex t — M^uJ^r 2 ^ + uj 2 z 2 )/2. In 
the limit of highly elongated traps (ui± ^> ui z ) the tight 
confinement ensures that no excited states are available 
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in the transverse direction and thus the dynamics takes 
place along the axial direction. Factorizing tp m (r) into a 
longitudinal and a transverse function, and approximat- 
ing the transverse part as the two-dimensional ground 
state of the transverse oscillator [ijl |25| , the equations 
of motion J3J become one-dimensional for the longitu- 
dinal wave functions i^ m (z), and the coupling constants 
c and c-i are accordingly rescaled by a factor l/(27rOj_), 
with a±_ = \J h/muj± the transverse oscillator length [30| • 

NUMERICAL METHOD 

The dynamics of the spin components is obtained by 
numerically integrating the coupled nonlinear differen- 
tial equations ©■ For the time evolution our numerical 
procedure combines the split operator method with the 
Fast Fourier Transform to treat the kinetic terms and 
a fourth-order Runge-Kutta method for the remaining 
terms of the dynamical equations. We have compared 
our combined numerical method with an evolution us- 
ing a pure fourth-order Runge-Kutta as the one used in 
Ref. [23j , and we have found that our method allows larger 
time steps making the computation notably faster. 

In this paper, we consider N = 20000 atoms of spin-1 
87 Rb trapped in a quasi-lD harmonic trap with oj± = 
2n x 891Hz and lu z = 2n x 21 Hz. The coupling con- 
stants are clq = 101. 8as, and — 100. 4as |3l| . with as 
the Bohr radius. Since C2 < the atomic interactions 
for 87 Rb atoms are ferromagnetic. The initial wave func- 
tions ip m (z,t = 0) have the same spatial profile as the 
ground state of the scalar condensate with coupling con- 
stant co/(27ra^). Each initial wave function component 
is however normalized to the corresponding initial popu- 
lation in spin-m, N m . Since i/j m (z,t) are complex func- 
tions, ip m (z,0) = \ij) m (z, 0)| exp (i9 m ), the initial phases 
have also to be fixed. However, for initial symmetric con- 
figurations (Ni — N-i) the dynamical evolution depends 
only on one initial relative phase 9 = 20q — 6\ — 9-\ [23| . 
and thus only one initial phase has to be fixed. Experi- 
mentally the initial phases are not well determined and 
the results are averaged over several repeated measure- 
ments \u\ . Therefore we will average the dynamics over 
different initial relative phases randomly distributed over 
(0, 2tt) in order to obtain results easy to compare with 
experiments. 

SPINOR DYNAMICS AT T = 

We consider that initially a quasi-pure condensate in 
the m = spin component is populated. Spin compo- 
nent mixing requires, at least, a small seed of atoms 
populating the other components but keeping the to- 
tal magnetization equal to zero. In Figure ^ we plot 
the population of each spin component as a function 



of time at zero temperature, for the initial populations 
(Afi/iV,7Vo/7V,iV_i/iV) = (0.5%, 99%, 0.5%). In absence 
of an external magnetic field gradient the dynamical evo- 
lutions of the spin m = ±1 components are symmetric. 
Dashed lines correspond to the dynamical evolution with 
a given initial relative phase 6 — 0, and solid lines are 
the average over 20 random initial relative phases. The 
oscillations of the populations with a given initial phase, 
are smeared out by the average over different runs. It 
is interesting to point out the damping of the oscilla- 
tions obtained in our numerical calculations which is a 
finite size effect related to collapse phenomena charac- 
teristic of discrete anharmonic spectra [^3, and to de- 
phasing of Josephson oscillations j^. During the time 
evolution the magnetization is conserved, as it has been 
experimentally observed 0|. The populations oscillate 
around the ground state configuration of the system with 
M. = that in absence of applied external magnetic 
field is (25%, 50%, 25%) 0. Our numerical results are 
in qualitative agreement with the experimental measure- 
ments of Ref. [l5j obtained in a strongly anisotropic disk- 
shaped spin-1 condensate, where the relaxation to the 
steady state is also not monotonic but damped. 

In Fig|3we plot the density profiles of the spin com- 
ponents for the configuration 8 = (see FigQ. At the 
initial stages of the evolution, the population of the spin- 
component decreases due to the spin exchange interac- 
tion, and thus the ±1 spin components start to be popu- 
lated by the same amount, keeping the symmetry of the 
initial state. The total magnetization is conserved along 
the time evolution. Initially (i < 100 ms), the conversion 
of atoms from to ± 1 states mainly occurs at the cen- 
tral part of the condensate, where the density is higher 
and thus the coupling between different spin components 
is more effective, see Eq. (@J. 

Then, the ± 1 spin components swing back to the 
component, and vice versa. The oscillations between 
the populations of the m = and m = ±1 states are 
not regular and present a dynamical instability around 
t ~ 100 ms, when the large amplitude oscillations be- 
come small amplitude oscillations [2 HE 0- At this 
moment the condensate starts the multi-domain forma- 
tion process into small dynamical spin domains. A sim- 
ple estimation of the time scale for the appearance of the 
instability tdom ~ 2irh/ fan) has been provided b y st udy- 
ing the normal excitation modes of the system |35|, |36( . 
In our case, taking n~2x 10 14 cm~ 3 at the center of 
the trap, tdom ~ 140 ms, which is in agreement with our 
results. Notice that, in the present study, the analysis is 
performed directly from the spatial and temporal evolu- 
tion of the population of the spin components. 

The transfer of population remains coherent along all 
the time evolution, and it does not become chaotic even 
for large times. The number of small spin domains does 
not grow indefinitely, but it is limited by a characteristic 
size of the spin domains (Idom) that depends on the in- 
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ternal coupling between different spin components, and 
more weakly on temperature. The formation of multi- 
ple spin domains, also found recently in Ref . [3^, |3(| , is a 
result of the ferromagnetic character of the spin-1 87 Rb 
condensate: C2 < favors the spatial separation of ±1 
atoms from 0-atoms, as obtained in our numerical results. 
Since in our simulations the initial condition for m = ±1 
are identical, the density profiles of m = ±1 components 
are equal along the time evolution. Apparently, the value 
of Idom ~ 5 ^tm, being much smaller than the harmonic 
oscillator length, does not depend on the external trap- 
ping potential, but it is an intrinsic characteristic of the 
spin coupling strength [3bT |. 

The dynamical evolution between spin components can 
be easily interpreted in terms of the effective potential 
V£ff(z, t) felt by each spin component ip m {z, t), and from 
the continuity equations. For the same initial configura- 
tion as in Fig|3 we plot in Fig. [3| the effective poten- 
tials and * = V_ x (top panels), and the local 
transfer of population 5ho(z, t) and 8hi(z,t) = 8h-x(z,t) 
(bottom panels) for t — 0, 40 and 80 ms. For a fixed time 
i, 8h m (z,t) represents the local exchange of atoms be- 
tween spin components. If 8n m {z, t) < at that position 
the population of the spin-m component is decreasing 
and thus the atoms convert to the other spin compo- 
nents. For example at the first stages of the time evo- 
lution (t < 60 ms) 5hi(z, t) = 8h-i(z,t) > 0, and con- 
sequently Sho(z,t) < 0: the spin-exchange interaction is 
favoring the conversion from to ±1 states as it is shown 
in Fig^ At t = this is also true, but it cannot be ap- 
preciated from the scale of the figure. Moreover, since 
the minimum and maximum of 8fio{z,t) and Sh±i(z,t), 
respectively are at the center of the condensate, the spin- 
exchange mainly occurs at the central region, as we have 
already commented in Fig[5] At t = 80 ms 8n±i is posi- 
tive at the boundaries of the condensate, and negative at 
the central region, whereas 8riQ has the opposite behav- 
ior. Therefore, the population with m — ±1 is decreasing 
at the center and increases at the boundaries, and the ±1 
atoms convert to state mainly at the center, where the 
m = 0-condensate develops a new central peak. 

We have also performed simulations starting from 
other initial conditions. In particular, we observe that 
multi-domain formation is inhibited if the starting con- 
figuration coincides with the ground state composition 
(25%, 50%, 25%). Moreover, we find a good agreement 
with the very recent experimental results of Chang et al. 
[8?j |. starting from (0, 3/4, 1/4) and converging to (1/5, 
2/5, 2/5). 

SPINOR DYNAMICS AT FINITE T 

The spinor dynamics of these multicomponent quan- 
tum gases is rather complex due to the internal coupling 
between different spin components. We consider now 



thermal effects in the spinor dynamics. At low temper- 
ature, thermal excitations can be described within the 
Bogoliubov-de Gennes theory. Recently, finite temper- 
ature effects in the equilibrium density distribution of 
the condensed and non-condensed components of spin- 
1 trapped atoms has been investigated within Hartree- 
Fock- Popov theory We will investigate thermal 

effects in the spinor dynamics using a Bogoliubov-de 
Gennes description of the thermal cloud. For a highly 
elongated trap, the condensate is quasi-lD, and total 
density fluctuations are strongly suppressed at small tem- 
peratures, whereas phase fluctuations are relevant and in 
the Thomas-Fermi regime can be described analytically 
in terms of Legendre polynomials Pj(z) [13, Anal- 
ogously, we assume that spin fluctuations can be disre- 
garded in a first approximation. 

In Refs. 0, H3 it has been shown that at finite temper- 
ature each spin component has its own thermal cloud. 
Thus, we assume that initially the condensate corre- 
sponding to each spin-m component can be described by 
an order parameter tp m (z) = y / n rn (z) exp(z0 m ), which 
has a random fluctuating phase :39| 

*-« = " [fa 

(^ARtf) 2 )] 1/2 Pj(z/R TF ) (af + af*) (6) 

where n [z) is the equilibrium total density profile of 
the initial condensate, [i and Rtf are its chemical po- 
tential and Thomas- Fermi radius, 6j = hw 7 .\/jlj + l)/2 
is the spectrum of low lying axial excitations |40j, and 
a T ( a T*) are com pl ex amplitudes that replace the quasi- 
particle annihilation (creation) operators in the mean- 
field approach. In the numerical calculation in order 
to reproduce the quantum statistical properties of the 
phase |26f . a™ and a™* are sampled as random vari- 
ables with a zero mean value and (|a™| 2 ) = iVj", where 
NJ 1 = [exp(ej /ksT) — 1] _1 is the occupation number for 
the quasi-particle mode j 0, 0] . 

In Fig^l we plot the dynamical evolution of the spin 
populations at T — 0.2T C for the same initial populations 
as in FigQ] The critical temperature of Bose-Einstein 
condensation, T c = fku z N/ ln(2A), corresponds here to 
a single component Bose gas in a harmonic trap of fre- 
quency lu z . Solid lines correspond to the numerical re- 
sults averaged over 20 random initial relative phases, and 
dashed lines to a single run with initial phase 9 = 0. 
The interaction of the condensate atoms with the ther- 
mal clouds smears out the oscillations present at T = 
(FigD) and leads to an asymptotic configuration with all 
m components equally populated (equipartition) . 

Thermal effects in the density profiles of the spin com- 
ponents are shown in Fig0 The occurrence of phase 
fluctuations due to thermal excitations at t = Q trans- 
forms to modulations of the density during the time evo- 
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lution as in a sp in-polarized single component elongated 
condensate |26j. In a spinor condensate this leads to 
multi-domain formation at much earlier times than at 
T = (see Fig|2J|, and the density profiles are no longer 
symmetric at finite T. Moreover, the existence of differ- 
ent thermal clouds for each spin component [2(| breaks 
also the symmetry of the m — ±1 spin components of 
the initial state of system. Therefore the m = ±1 den- 
sity profiles are different during the time evolution, and 
the three components separate in different spin domains. 
The local magnetization is not longer conserved as at 
T = but, as expected, the total magnetization is still a 
conserved quantity along all the time evolution. Similar 
results are obtained at lower temperatures T/T c = 0.01 
and 0.001. 



CONCLUSIONS 



domain formation go along with a loss of phase relations, 
and give rise to some enhanced (generalized) phase fluc- 
tuations? These questions go beyond the present study, 
and will be discussed in detail in a future publication, in 
which we will compare the F = 1 and F = 2 case. 
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In this paper, we have studied the spinor dynamics of 
a spin-1 87 Rb condensate in a highly elongated trap. We 
have solved the full three coupled dynamical equations 
for the spin components within Gross-Pitaevskii frame- 
work without any further approximations. This is in 
fact necessary, since approximated approaches frequently 
mask some of the aspects of the dynamics. We have also 
considered here finite temperature effects, using the ap- 
proach of Ref . 01 ■ We have found that the spinor dy- 
namics towards the steady state is not monotonic but 
rather slowly damped, involving a coherent transfer of 
population between different spin components. At finite 
temperature the coherent oscillations of populations are 
almost smeared out. The internal coupling of the spin 
components leads to the formation of multiple spin do- 
mains of a small, but finite characteristic length Idom 
which does not decrease with time, and seems thus to 
be determined intrinsically by the nonlinear interactions. 
This scale is evidently larger than the condensate heal- 
ing length, Iheal = 2nh/^2Mc n, which for scattering 
lengths a of order of 5 nm is of order of 10-100 nm. 
In fact, Idom is of the order of the spin healing length 
£ s = 2-nh/ W2M\ca\n. The presence of different thermal 
clouds for each spin component breaks the symmetry of 
the m = ±1 components, and therefore separates them 
in different spin domains. For a condensate with initially 
zero magnetization, the spin populations oscillate around 
the ground state configuration (25%, 50%, 25%) at T = 0, 
whereas at finite temperature the interaction of the con- 
densate atoms with their corresponding thermal clouds 
leads to equipartition in populations, i.e. (1/3, 1/3, 1/3). 

Our results shed, in our opinion, also some light on the 
question of decoherence. Our simulations, and in partic- 
ular the finding of multi-domain formation suggest that 
decoherence undergoes enhancement with the number of 
components in the system. Of course, there are many 
open questions connected to this, e.g. does the multi- 
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FIG. 1: (Color online) Population of the spin compo- 
nents as a function of time for the initial configuration 
(Ni/N,N /N,N-!/N) = (0.5%, 99%, 0.5%) at T = 0. Solid 
lines: numerical results averaged over 20 random initial rela- 
tive phases 9. Dashed lines: initial phase 9 = 0. 
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FIG. 2: (Color online) Density profiles of the spin compo- 
nents at different times (in ms) at T = 0, ni(z,t) (black) 
and no(z,t) (green). The initial configuration corresponds to 
(0.5%, 99%, 0.5%), and 9 = (dashed line in FigQ. 
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FIG. 3: (Color online) For the initial configuration 
(0.5%, 99%, 0.5%), with 6 = and T = 0. Top panels: Ef- 
fective potentials of the spin components for t = 0, 40 and 80 
ms. Bottom panels: Transfer of populations 8n m (z,t) for the 
same times. 



0.75 



0.25- 







i 

T= 0.2 T 


\\ 






V- 




/ ' *~ 


_y i 




1,1, 



100 200 300 400 500 

t (ms) 



FIG. 4: (Color online) Population of the spin components as a 
function of time for the initial configuration (0.5%, 99%, 0.5%) 
at T = 0.2T C . Solid lines: numerical results averaged over 20 
random initial relative phases 0. Dashed lines: initial phase 
9 = 0. 
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FIG. 5: (Color online) Density profiles of the spin components 
at different times (in ms) at T — 0.2T C for the same initial 
configuration as in Figl2l 



